Trait Ordination Plant Community Analysis

Trait Ordination Plant Community Analysis

1. traits

   1. above- and belowground

      1. height & SLA
         

   2. acquisitive vs. conservative resource-use strategies

      1. rooting depth

   3. collaborative vs. individual resource-use strategies

      1. fine root proportion

      2. thickness of roots

Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.

Author: Caryn D. Iwanaga

Updated: 01/24/2025

Libraries

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-8

Read in Data

cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
  mutate(
    nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
    ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))


# Define a reusable color scale
# custom_colors <- scale_fill_manual(
#   values = c(
#     "control" = rgb(195, 197, 193, maxColorValue = 250),
#     "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
#     "+compost" = rgb(167, 156, 109, maxColorValue = 225)
#   )
# )

custom_colors <- scale_fill_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

# Set a global theme
theme_set(theme_bw())

Make site x species matrix

1. Calculate relative abundances of species for each site
   1. confused how to aggregate abundances because even grouping by plotid (nut_trt and ppt_trt), year, species, and grazing history (block 1&2 – high, block 3&4 – low) and averaged, the sd are as high as 45.96!
      1. proceeded with these high sd to see what happens (ask Lauren)

site.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  mutate(
    grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
  ) %>% # separate by block (grazing intensities)
  group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

matrix0 <- spread(
  site.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix <- matrix0[, c(5:79)]

# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)

Run NMDS

Rules of thumb:

• < 0.05 is excellent

• 0.05-0.1 is good

• 0.1-0.2 is fair

• 0.2-0.3 is cause for concern…

• > 0.3 is poor, random

All years

Stress: fair

0.1966553 

Grouping by year and grazing history had most significant clustering shown.

metaMDS()

nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)

Run 0 stress 0.2020736 
Run 1 stress 0.2023246 
... Procrustes: rmse 0.07023633  max resid 0.190107 
Run 2 stress 0.2017185 
... New best solution
... Procrustes: rmse 0.02245434  max resid 0.1414201 
Run 3 stress 0.2088616 
Run 4 stress 0.1970656 
... New best solution
... Procrustes: rmse 0.07244987  max resid 0.2329502 
Run 5 stress 0.2176689 
Run 6 stress 0.1970626 
... New best solution
... Procrustes: rmse 0.006282932  max resid 0.03757378 
Run 7 stress 0.2020248 
Run 8 stress 0.2200448 
Run 9 stress 0.2033604 
Run 10 stress 0.1973254 
... Procrustes: rmse 0.02289426  max resid 0.1158253 
Run 11 stress 0.2042184 
Run 12 stress 0.2020728 
Run 13 stress 0.2017572 
Run 14 stress 0.203447 
Run 15 stress 0.2183276 
Run 16 stress 0.1990614 
Run 17 stress 0.2000129 
Run 18 stress 0.2089696 
Run 19 stress 0.2082397 
Run 20 stress 0.2028091 
Run 21 stress 0.2000005 
Run 22 stress 0.1999666 
Run 23 stress 0.2019172 
Run 24 stress 0.2147543 
Run 25 stress 0.2017483 
*** Best solution was not repeated -- monoMDS stopping criteria:
    25: stress ratio > sratmax

nmds.comp # pretty not great stress

Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.1970626 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 6 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray’ 

stressplot(nmds.comp)

plot(nmds.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))

nmds_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years - LOW GRAZING ONLY

Stress: fair

0.1869068

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
  filter(grazing_hist == "low")

row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")

matrix.low <- matrix.low[, c(5:79)]

rownames(matrix.low) <- row.low

Warning: Setting row names on a tibble is deprecated.

matrix.bray.low <- decostand(matrix.low, "total")

metaMDS()

plot

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  color = as.factor(yr), shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

# add year ellipses ---- 
ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
    stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

Make each individual year NMDS to see if there’s better stress

2019

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
  filter(yr == 2019)

row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")

matrix.2019 <- matrix.2019[, c(5:79)]

rownames(matrix.2019) <- row.2019

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2019 <- decostand(matrix.2019, "total")

nmds 2019

Stress: Excellent

0.0560178

metaMDS()

nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)

Run 0 stress 0.0560178 
Run 1 stress 0.05957995 
Run 2 stress 0.05960972 
Run 3 stress 0.05851815 
Run 4 stress 0.05851812 
Run 5 stress 0.05957995 
Run 6 stress 0.06407773 
Run 7 stress 0.05851805 
Run 8 stress 0.05982901 
Run 9 stress 0.06074817 
Run 10 stress 0.05984339 
Run 11 stress 0.05982901 
Run 12 stress 0.07612032 
Run 13 stress 0.06031567 
Run 14 stress 0.05652852 
Run 15 stress 0.06130516 
Run 16 stress 0.05905532 
Run 17 stress 0.0560178 
... Procrustes: rmse 1.328847e-06  max resid 3.511506e-06 
... Similar to previous best
Run 18 stress 0.07356079 
Run 19 stress 0.05652852 
Run 20 stress 0.05957995 
*** Best solution repeated 1 times

nmds.2019 # pretty not great stress

Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2019 
Distance: bray 

Dimensions: 2 
Stress:     0.0560178 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2019’ 

stressplot(nmds.2019)

plot(nmds.2019, type="t")

nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
  )
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))

nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

# read in trait data

# subset trait data to match comp data

# set matching row names

2020

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
  filter(yr == 2020)

row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")

matrix.2020 <- matrix.2020[, c(5:79)]

rownames(matrix.2020) <- row.2020

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2020 <- decostand(matrix.2020, "total")

nmds 2020

Stress: fair

0.1522495

metaMDS()

nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)

Run 0 stress 0.1524241 
Run 1 stress 0.1729906 
Run 2 stress 0.1595535 
Run 3 stress 0.1584809 
Run 4 stress 0.1524241 
... New best solution
... Procrustes: rmse 1.001013e-05  max resid 2.563444e-05 
... Similar to previous best
Run 5 stress 0.1522495 
... New best solution
... Procrustes: rmse 0.02107748  max resid 0.06346541 
Run 6 stress 0.1524241 
... Procrustes: rmse 0.02107671  max resid 0.06351798 
Run 7 stress 0.1522495 
... New best solution
... Procrustes: rmse 6.786807e-06  max resid 1.943279e-05 
... Similar to previous best
Run 8 stress 0.1821737 
Run 9 stress 0.1594876 
Run 10 stress 0.1872714 
Run 11 stress 0.1705568 
Run 12 stress 0.1752228 
Run 13 stress 0.1522045 
... New best solution
... Procrustes: rmse 0.01182953  max resid 0.03730308 
Run 14 stress 0.1747777 
Run 15 stress 0.1592892 
Run 16 stress 0.1595535 
Run 17 stress 0.1625903 
Run 18 stress 0.1594876 
Run 19 stress 0.165848 
Run 20 stress 0.1594877 
Run 21 stress 0.1524241 
... Procrustes: rmse 0.0248137  max resid 0.06599096 
Run 22 stress 0.1584809 
Run 23 stress 0.1522495 
... Procrustes: rmse 0.01183773  max resid 0.03738525 
Run 24 stress 0.1744238 
Run 25 stress 0.1744238 
*** Best solution was not repeated -- monoMDS stopping criteria:
    16: stress ratio > sratmax
     9: scale factor of the gradient < sfgrmin

nmds.2020

Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2020 
Distance: bray 

Dimensions: 2 
Stress:     0.1522045 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2020’ 

stressplot(nmds.2020)

plot(nmds.2020, type="t")

nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
  )
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))

nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 1) + 
  labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment") 

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Nutrient Treatment")

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Precipitation Treatment")

2021

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
  filter(yr == 2021)

row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")

matrix.2021 <- matrix.2021[, c(5:79)]

rownames(matrix.2021) <- row.2021

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021 <- decostand(matrix.2021, "total")

nmds 2021

Stress: fair

0.1230831

metaMDS()

nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)

Run 0 stress 0.1232098 
Run 1 stress 0.132972 
Run 2 stress 0.132972 
Run 3 stress 0.1506988 
Run 4 stress 0.1515614 
Run 5 stress 0.1232098 
... New best solution
... Procrustes: rmse 1.844754e-06  max resid 4.667063e-06 
... Similar to previous best
Run 6 stress 0.1230832 
... New best solution
... Procrustes: rmse 0.04464562  max resid 0.1321501 
Run 7 stress 0.1232098 
... Procrustes: rmse 0.04464457  max resid 0.1315388 
Run 8 stress 0.1421227 
Run 9 stress 0.1339879 
Run 10 stress 0.1515617 
Run 11 stress 0.1421227 
Run 12 stress 0.1506988 
Run 13 stress 0.132972 
Run 14 stress 0.1555173 
Run 15 stress 0.1436877 
Run 16 stress 0.1761755 
Run 17 stress 0.1339879 
Run 18 stress 0.1232098 
... Procrustes: rmse 0.04464533  max resid 0.1315293 
Run 19 stress 0.1507047 
Run 20 stress 0.1525517 
Run 21 stress 0.1339879 
Run 22 stress 0.1230832 
... Procrustes: rmse 2.173671e-05  max resid 6.138502e-05 
... Similar to previous best
*** Best solution repeated 1 times

nmds.2021

Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021 
Distance: bray 

Dimensions: 2 
Stress:     0.1230832 
Stress type 1, weak ties
Best solution was repeated 1 time in 22 tries
The best solution was from try 6 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021’ 

stressplot(nmds.2021)

plot(nmds.2021, type="t")

nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
  )
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))

nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Precipitation Treatment")

Make species x trait matrix

Greenhouse traits

• height, leaf dry matter content, specific leaf area, root mass fraction, relative growth rate, coarse root diameter, root density, specific root length of fine roots, specific root length of coarse roots

Field traits

• height, leaf dry matter content, specific leaf area, seed mass, leaf C, leaf N

# read in trait data

# subset trait data to match comp data

# set matching row names